326 8.2  Molecular Simulation Methods

Similarly, the occupancy probability in each microstate can be estimated using a weighted

histogram analysis method, which enables the differences in thermodynamic parameters for

each state to be estimated, for example, to explore the variation in free energy differences

between the intermediate states.

An extension to Monte Carlo MD simulations is replica exchange MCMC sampling, also

known as parallel tempering. In replica exchange, multiple Monte Carlo simulations are

performed that are randomly initialized as normal, but at different temperatures. However,

using the Metropolis criterion, the molecular configurations at different temperatures

can be dynamically exchanged. This results in efficient sampling of high and low energy

conformations and can result in improved estimates to thermodynamical properties by

spanning a range of simulation temperatures in a computationally efficient way.

An alternative method to umbrella sampling is the simulated annealing algorithm, which

can be applied to Monte Carlo as well as other molecular simulations methods. Here the

simulation is equilibrated initially at high temperature but then the temperature is gradually

reduced. The slow rate of temperature drop gives a greater time for the simulation to explore

molecular conformations with free energy levels lower than local energy minima found by

energy minimization, which may not represent a global minimum. In other words, it reduces

the likelihood for a simulation becoming locked in a local energy minimum.

8.2.4  AB INITIO MD SIMULATIONS

Small errors in approximations to the exact potential energy function can result in artifac­

tual emergent properties after long simulation times. An obvious solution is therefore to

use better approximations, for example, to model the actual QM atomic orbital effects, but

the caveat is an increase in associated computation time. Ab initio MD simulations use the

summed QM interatomic electronic potentials UQD experienced by each atom in the system,

denoted by a total wave function ψ. Here, UQD is given by the action of Hamiltonian operator

Ĥ on ψ (the kinetic energy component of the Hamiltonian is normally independent of the

atomic coordinates) and the force FQD by the grad of UQD such that

(8.19)

U

r

H

F

U

r

QD

QD

QD

→

→

^

( ) =

−∇

( )

(

)

=

ψ

These are, in the simplest cases, again considered as pairwise interactions between the highest

energy atomic orbital of each atom. However, since each electronic atomic orbital is identified

by a unique value of three quantum numbers, n (the principal quantum number), ℓ (the azi­

muthal quantum number), and mℓ (the magnetic quantum number), the computation time for

simulations scales as ~O(n³). However, Hartree–​Fock (HF) approximations are normally applied

that reduce this scaling to more like ~O(n^2.7). Even so, the computational expense usually restricts

most simulations to ~10–​100 atom systems with simulation times of ~10–​100 ps.

The HF method, also referred to as the self-​consistent field method, uses an approximation of

the Schrödinger equation that assumes that each subatomic particle, fermions in the case of elec­

tronic atomic orbitals, is subjected to the mean field created by all other particles in the system.

Here, the n-​body fermionic wave function solution of the Schrödinger equation is approximated

by the determinant of the n × n spin–​orbit matrix called the “Slater determinant.” The HF elec­

tronic wave function approximation ψ for an n-​particle system thus states

(8.20)

ψ

χ

χ

χ

χ

χ

χ

→→

→

→

→

→

→

→

…

…

r r

r

n

r

r

r

r

r

n

n

n

1

2

1

1

1

1

1

1

2

2

2

1

, ,

,

…

(

) ^≈

( )

( )

( )

( )

( )

→

→

→

→

⋯

⋯

⋱

⋯

…

r

r

r

r

n

n

n

2

1

2

2

( )

( )

( )

( )

χ

χ

χ